Clustering Kinetics of Granular Media in Three Dimensions 
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Three-dimensional molecular dynamics simulations of dissipative particles (~ 10 6 ) are carried 
out for studying the clustering kinetics of granular media during cooling. The inter-connected high 
particle density regions are identified, showing tube-like structures. The energy decay rates as 
functions of the particle density and the restitution coefficient are obtained. It is found that the 
probability density function of the particle density obeys an exponential distribution at late stages. 
Both the fluctuation of density and the mean cluster size of the particle density have power law 
relations against time during the inelastic coalescing process. 
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To understand the fundamental differences of macro- 
scopic behaviors between a regular gas and the inherently 
dissipative granular media is a challenging problem. Re- 
cent studies 0] have revealed that inelastic collision in 
granular media gives rise to a variety of complex phe- 
nomena which are absent in classical ideal gas, including 
the formation of clustering structures |^-|5| and collaps- 
ing (^,0]. To understand the clustering kinetics and to 
make the connection between the microscopic and macro- 
scopic properties of granular media are not only crucial 
for constructing the constitutive equations of granular 
materials, but also important for various engineering ap- 
plications, including packing ||, size segregation PJToJ 
and transport of granular materials p]Jll[|. 

Cluster formation has been studied before. In partic- 
ular, using a model of hard disks with dissipation, Gold- 
hirsch and Zanetti found that the clustering parti- 
cles form long string- like structure in two dimensions. 
Experimentally, two dimensional cluster formation has 
been studied in a driven system by Kudrolli et al. 
using a system consisting of small steel spheres rolling 
on a smooth bounded rectangular surface with moving 
side walls. They measured the probability density func- 
tions (PDF) of the particle density and the particle ve- 
locity for the clusters near the driven wall. However, 
most existing numerical simulations and experiments for 
clustering dynamics are carried out for two dimensions 
and only a small number of granular particles (in the 
order of 10,000) are used. For such a system, it is dim- 
cult to obtain reliable statistics and to address problems 
connecting macroscopic description with microscopic dy- 
namics. 

In this Letter, we present an analysis of three- 
dimensional clustering kinetics of granular media in 
a free cooling condition using a recently developed 
three-dimensional molecular dynamics (MD) code based 
on Massage-Passing- Interface (MPI). This parallel code 
makes it possible to simulate over one million dissipative 
particles in three dimensions and therefore allows us to 
tackle problems which are difficult in previous studies, in- 



cluding the growth of fluctuation density and the mean 
cluster size of the particle clusters. In this paper, we 
focus on studying the energy decaying dynamics, three- 
dimensional structures of clusters and the scaling kinetics 
of the particle density. 

The granular media in this paper is treated as a collec- 
tion of identical spherical particles with radius R. The 
motion of each particle is simulated by Newton's second 
law while ignoring the angular motion. To specify the 
interaction forces in our MD simulation model, we as- 
sume that there are no interactions when particles are 
not in contact. The following interaction forces are in- 
cluded when two particles, i and j overlap |l2|, i.e., the 
distance |ry| (r.y = Tj — r^) is smaller than 2R: (1) An 

elastic restoration force: fj£ — ym,- (|ry| — 2iJ)rjj/|ry |; 
(2) dissipations due to the inelasticity of the collision: 
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Here Y is 
is the mass of particle i; 



v™- are the prO- 



^ar ~ 'Js™^^. 

the Young's modulus, rrii ~ R 3 
v"j = (vij ■ ry)rij/|ry| 2 and vf^ = v y - 
jections of the relative velocities (= v, — v?) on the 
direction and the tangential direction, respectively; 7 and 
7 S are the dissipation coefficients for the relative motion 
of particles in the normal and tangential directions. 

Initially particles are distributed uniformly in a cu- 
bic box with a small random perturbation. The initial 
particle velocities are drawn from an isotropic Gaussian 
distribution with variance a. For simplicity, in this paper 
we only report results for a = 4 and the Young's mod- 
ulus Y = 10 5 . The radius of granular particles, R, and 
the mass of the particles, to, are normalized to be 1. A 
second-order scheme is used for the time integration of 
the motion of particles and simulations are carried out 
in a periodic box of size L 3 . Because of the large value 
of Young's modulus Y, our time step At has to be set 
to At = 4.0 x 10 -4 , small enough to resolve the full col- 
lision process. The particle volume fraction a is defined 
as a — 4ttN/3L 3 , where N is the particle number in 
the system. The restitution coefficient, r, defined as the 
ratio of the relative speeds after collision and before col- 
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lision, can be determined as: r = \J (rj* + 2r 2 )/3, where 



r„ = exp(— tt/ -v/ oj 2 /t) 2 
and 7/ = j /2m J13|. 
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FIG. 1. A typical particle configuration at late stages of clus- 
tering formation for parameters a = 0.074 and r = 0.6. The 
simulation is for 10 particles, but only one-sixty-fourth of 
the simulation domain is shown. 




FIG. 2. Iso-surface of the particle density for p = 3.5p m . The 
resolution in this plot is 50 3 and the simulation was carried 
out for L = 384. 

In Fig. 1, we show a snapshot of a three dimensional 
particle distribution in a cooling granular system at late 
stage t = 1000 from a simulation of one million particles 
(the initial characteristic collision time is of 0(1)). The 
parameters in this calculation arc a = 0.074, 7 = 102 and 
7 S = 51, which result in r — 0.6. For clarity, only one- 
sixty-fourth of the whole simulation domain is shown in 
this plot. To calculate the local particle density /a(x), we 
equally divide the whole space into small sub-boxes and 
average the particle number in each box. In Fig. 2, we 
present iso-surfaces of the particle density function p(x) 
for p = 3.5p m , where p m is the mean density of the whole 



space. From Fig. 1 and Fig. 2, it is evident that the par- 
ticle clusters are formed and are more or less spatially 
inter-connected. In addition, the cluster structures are 
tube-like, which is to be compared with the string-like 
particle structure in two dimensions |^|. We have also 
calculated fractal dimensions for different regions using 
the box-counting method |l4| . We found that the fractal 
dimension df ~ 1.3 for p > 3.5p m , in a qualitative agree- 
ment with the observed tube- like structure in Fig. 2, and 
df ~ 2.2 for p > p m . 

In a cooling system, the forms of the particle motion 
fall roughly into two categories: free streaming and colli- 
sion. During the streaming process, the speed of each 
particle is constant. The collision step decreases the 
speed of particles involved in the collision due to dis- 
sipation and reduces the energy. Based on a kinematic 
argument, it has been shown in [§,|l5| that if the spatial 
particle distribution is uniform for which the mean free 
path is proportional to E 1 / 2 , and the local velocity dis- 
tribution is random which leads to strong collision with 
loss of energy proportional to E, then the total energy of 
the system, E(t) = 1/2^ v?, decays as: 



E(t) 



E(0) 



(1) 



(1 + t/t 

where t c is a crossover time, and the exponent (3 equals 
2. 

Our simulations show that in general the above equa- 
tion can be used to describe the decay of the total energy 
during the particle clustering process, but depends on 
the properties of the granular media, in particular, the 
particle volume fraction and the restitution coefficient. 
Fig. 3a shows the energy decay exponent (3 as a function 
of the particle volume rate a for the restitution coeffi- 
cient r = 0.43. In the inset, the total energy versus time 
for a typical case is shown. One can see that a power 
law decay can be identified for the late stage. Similar 
results are shown in Fig. 3b for against the restitution 
coefficient r with a = 0.09. The results shown here give 
a quantitative description of the dependence of the decay 
exponent on various material parameters: in the case of 
low density or large restitution, the particle distribution 
in the system is almost uniform and the resulting decay 
exponents equal 2, in a good agreement with the heuris- 
tic argument in [^|,^5). The smaller decay exponents are 
seen for systems with high densities and small restitu- 
tions, where clustering occurs. In this circumstance, the 
speed of particles in a cluster is approximately the same 
and the collision frequency is relatively low, leading to 
less energy dissipation and smaller energy decay rate. 

The decrease of the energy decay rate has been studied 
previously in using a linear hydrodynamic analysis 
based on phenomenological macroscopic equations. It is 
argued that for a system with particle clusters, the hy- 
drodynamic modes have low decay rates. Quantitatively, 
it is interesting to notice that the energy decay exponent 
[3 is round 1.4 ~ 1.5 in the high density and small resti- 
tution regimes. To understand this nontrivial exponent, 
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we decompose the total energy: E(t) = E^(t) + Ef(t), 
with Ek(t) and Ef(t) being the local spatially averaged 
kinetic energy and the fluctuation energy, respectively. 
We argue that in the high density or small restitution 
case, due to the formation of particle clusters, the mo- 
tion of the clustered particles dominates, therefore the 
decay rate of E(t) follows directly the macroscopic hy- 
drodynamic decay of Ek(t): 



E k {t)~ J \v(k,t = 0)\ 2 exp(-2v\k\ 2 t)d D k 



f 



-D/2 



(2) 



where v is some effective viscosity and D is the spatial 
dimension. In deriving (2), We have assumed that the 
initial velocity fluctuation v(k, t = 0) is independent of 
k for small \k\ (certainly true for the initial velocity dis- 
tribution). For D = 3, Our numerical results of f3 ~ 1.5 
in the clustering regime agree well with Eq. (2). To fur- 
ther verify our prediction, we have also carried out 2D 
simulations (not shown here). We found that the energy 
decay rate at high densities is indeed close to 1, while the 
exponent /3 for very low density remains 2. 
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FIG. 3. The energy decay exponent (5 as functions of the 

particle volume fraction (a) and the restitution coefficient (b). 

The inserted pictures show the typical energy decays against 

time for (a) a = 0.011; (b) r = 0.20. 

In the absence of clustering, the spatial particle distri- 
bution should be uniform and the PDF of the particle 
density is given by a <5-function. On the other hand, 
if clustering occurs, particles coalesce to certain regions 
while leaving others void, leading to the disparity of par- 
ticle density distribution. The study of the time depen- 
dence of the PDF of the particle density greatly enhances 



our knowledges of the clustering process and the charac- 
teristics of the clusters. It is interesting to note that the 
shape of the PDF is correlated with the energy decay 
exponent (3. For j3 ~ 2, the corresponding PDF of the 
density shows very little spreading, signaling the non- 
existence of clustering in that parameter range, whereas 
in the clustering parameter regime (e.g., the parameters 
used in Fig. 4), j3 is close to 1.5. In Fig. 4, we present the 
PDF of the density function p at different times with the 
same parameters as in Fig. 1. The statistics are obtained 
using 32 3 sub-boxes. The density is normalized by the 
average density p m . The PDF at t = is centered at 
one narrow density interval due to our specific initializa- 
tion of the particle distribution. As the systems evolves, 
the PDF of the density function widens, implying the 
non-uniformity of the spatial particle distribution and 
the formation of clusters. Finally, the system becomes 
totally spread with more samples at lower density and 
fewer at high density. The PDFs at late stages can be 
approximated by an exponential function, indicating the 
existence of all possible densities. No local peak of PDF 
at high density was observed, at least for the simulation 
times we have reached, implying no specific characteris- 
tic density scale. Finally, in the inset of Fig. 4, we show 
the density fluctuation 5p(t) as a function of time, which 
follows a power law for almost three decades in time: 



0.55 



(3) 
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Normalized Density p/p m 
FIG. 4. Normalized PDFs of the particle density at differ- 
ent times. The inserted picture shows that the fluctuation of 
the density function grows as a power function of time. 

To study the particle density correlation, we ana- 
lyze the power spectra of the particle density function: 

E p( K ) = T,K-l/2<k<K+l/2p( k )P*( k )> whcrc P( k ) is 

the particle density in Fourier space. In the inset of 
Fig. 5, E p {K) as a function of K at different times are 
shown. Two features in the power spectra can read- 
ily be seen. First, as the time evolves, the peak value 
of the power spectra increases, in agreement with the 
inset in Fig. 4 due to the fact 8p 2 = J K E p (K)dK. 
Second, if we denote K max (t) for the wavenumber at 
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which E p {K) has the maximum value, we can see that 
K max (t) decreases as a function of time, indicating the 
increase of the typical cluster size. To quantify this be- 
havior, we define a mean cluster size, X(t) by: X(t) = 
J K E p (K)dK/ f K E p (K)KdK. In Fig. 5, we present \{t) 
as a function of time for the times when particles form 
clusters. Again, a power law dependence of X(t) on t is 
observed over two decades in time [fill : 

A(t)~t< , (4) 

where £ is the growth exponent for the mean cluster size 
and is close to 0.13 for the given case. From other simu- 
lations not shown here, we note that C has a week depen- 
dent on the particle volume fraction and increases slowly 
with increasing of a. 
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FIG. 5. The mean size of particle cluster as a function of 
time. In the inset shows the power spectra E P (K) at different 
times. 

To summarize, we have studied the clustering phe- 
nomenon in the cooling process of three dimensional 
granular media. Large spatial density fluctuation is ob- 
served due to inelastic collision and the formation of the 
clusters strongly affects the energy decay rate. It is found 
that the total energy has a power law decay against time 
at late stages: E(t) ~ t~P with [3 ~ 1.5 in the parameter 
regimes where clustering forms. The value of (3 crosses 
over to its mean field value (3 m ft = 2.0 for systems with 
large restitution constant r or small density a where the 
particle clustering is absent. By making simple analogy 
to the decay of kinetic energy in fluids, we relate the 
value of (3 in the clustering parameter regime to the spa- 
tial dimension D: (3 = D/2. We observe that clusters in 
three dimensions form tube-like structures for high par- 
ticle density regions at late stages. A preliminary study 
has revealed that the structure might be fractal. We 
have also studied the dynamical properties of clusters. 
It is found that the density fluctuation and the cluster 
size both obey power law scalings against time with non- 
trivial scaling exponents. We have no theory to explain 



the observed scaling dynamics and it is our hope that the 
present results would stimulate further studies in under- 
standing these intriguing behaviors. 
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